Import data

In this document we are going to calculate Evasion from the Vaisala Sensors and compare them to the eosFD flux measurements as well as see how they relate to other variables such as precipitation and temperature. First we import the data and for now we will filter so we are looking at non-injections

Create new time variables for ease of plotting

These new variables will help us plot each seperate day on a timescale of midnight -> midnight

df$day <- as.factor(substr(as.character(df$DateTime),1,10))
times <- substr(as.character(df$DateTime),12,16)
df$time <- strftime(df$DateTime, format="%H:%M")
df$time <- as.POSIXct(df$time, format="%H:%M")

Change in CO2 Station 1->4

Let’s visualize the changes in diel CO2 concentrations between a pair of sensors. Here we look at the entire reach, so stations 1 -> 4. We have a filter to remove some data that had issues.

filt <- df%>%
  filter(!is.na(V1)&!is.na(V4))%>%
  filter(DateTime > as.POSIXct("2019-07-13 00:00:00"))%>%
  filter(DateTime < as.POSIXct("2019-08-13 19:00:00")|DateTime > as.POSIXct("2019-08-14 23:59"))%>%
  filter(DateTime < as.POSIXct("2019-07-14 00:00")|DateTime > as.POSIXct("2019-07-14 23:59"))%>%
  filter(!DateTime == as.POSIXct("2019-08-09 09:57:00"))


dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(-3200,-3200,-100,-100))

my_breaks <- c(0,15,50,166) # Create breaks

plot14Ppm <- ggplot(filt)+
  geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
  geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
  geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
  geom_point(aes(x = time, y = V4-V1,group = day, col=stn1_Q),lwd=2)+
  scale_color_gradient(name = "Discharge (L/s)", trans = "log",
                        breaks = my_breaks, labels = my_breaks)+
  labs(x = "Time",y="DCO2", title = "DCO2 Change station 4 minus station 1")+
  theme_bw()

plot14Ppm

So now that we know what the change in CO2 between stations, it’s time to account for stream metabolism. Metabolism can be determined by taking the change in CO2 between stations before sunrise (photosynthesis is zero) and subtracting it from the change in CO2 at all other times. Let’s try to confirm this by looking the changes in DO over the same period.

Photosynthesis: Use CO2 and create O2 (1:1 ratio).

Respiration: Use O2, create CO2 (1:1 ratio).

Respiration can occur in the abscence of sunlight while photosynthesis cannot. At night, stream processes are respiration driven. During the day both respiration and photosynthesis occur. Is respiration constant throughout the day? if it is then pre-sunlight change in DO can be sen as respiration if no evasion of O2 is occuring. Therefore the area under the curve of the change in CO2 would be photosynthesis - respiration.

Metabolism

Metabolism

Change in Dissolved Oxygen station 1 -> 2

dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(-1,-1,1,1))

plot12DO <- ggplot(filt)+
  geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
  geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
  geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
  geom_point(aes(x = time, y = DO2_mg.L-DO1_mg.L,group = day, col=tempC_436),lwd=2)+
  scale_color_gradient(name = "Discharge (L/s)", trans = "log",
                        breaks = my_breaks, labels = my_breaks)+
  labs(x = "Time",y="DO (mg/L)", title = "Station 2 DO minus Station 1 DO")+
  theme_bw()

plot12DO
## Warning: Removed 3136 rows containing missing values (geom_point).

Change in Dissolved Oxygen station 1 -> 4

dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(0,0,4.5,4.5))

plot14DO <- ggplot(filt)+
  geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
  geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
  geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
  geom_point(aes(x = time, y = DO4_mg.L-DO1_mg.L,group = day, col=stn1_Q),lwd=2)+
  scale_color_gradient(name = "Discharge (L/s)", trans = "log",
                        breaks = my_breaks, labels = my_breaks)+
  labs(x = "Time",y="DO (mg/L)", title = "Station 4 DO minus Station 1 DO")+
  theme_bw()

plot14DO
## Warning: Removed 3136 rows containing missing values (geom_point).

DO vs CO2

filt%>%
  plot_ly(x = ~V1, y = ~DO1_mg.L,hoverinfo = "text",
        text = ~paste("Time: ", DateTime,
        "<br> Discharge: ", round(stn1_Q,2), "liters per second"))%>%
  add_markers(color = ~tempC_436)%>%
  layout(
    title = "Station 1 Dissolved Oxygen vs CO2",
    xaxis = list(title = "CO2 (ppm)", type = "log"),
    yaxis = list(title = "Dissolved Oxygen (mg/L)")
  )
## Warning: Ignoring 3136 observations

DO vs Temperature

## Warning: Ignoring 3345 observations

DO at each sensor

## Warning: Removed 3136 rows containing missing values (geom_point).

## Warning: Removed 3136 rows containing missing values (geom_point).

What’s up with station 4 DO

There is a diel drop in DO in the late afternoon that seems to be related to temperature

## Warning: Removed 3136 rows containing missing values (geom_point).

Dissolved Oxygen and Water Temperature

(Station 4)

## Warning: Ignoring 3345 observations

Calculate respiration between stations 1 & 4

Find the max change in pCO2 which will be when respiration is the least (before sunrise). This level serves as the baseline for that day (zero uptake). For every other time, we take the value of V1-V4 and subtract the zero uptake value for that day from it to determine uptake.

filtsub <- filt%>%
  filter(time < as.POSIXct("06:30", format = "%H:%M"))%>%
  group_by(day)%>%
  mutate(baseDeltCo2 = max(V1-V4))%>%
  ungroup()%>%
  select(day,baseDeltCo2)%>%
  distinct()

filt <- left_join(filt,filtsub)%>%
  mutate(Uptake = baseDeltCo2-(V1-V4))
## Joining, by = "day"
dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(-500,-500,1800,1800))

plotUptake <- ggplot(filt)+
  geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
  geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
  geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
  geom_point(aes(x = time, y = Uptake,group = day, col=stn3_Q),lwd=2)+
  scale_color_gradient(name = "Discharge (L/s)", trans = "log",
                        breaks = my_breaks, labels = my_breaks)+
  labs(x = "Time",y="Uptake", title = "CO2 Respiration station 1 -> 4")+
  theme_bw()
plotUptake
## Warning: Removed 78 rows containing missing values (geom_point).

Calculate Evasion using my equation

Now that we know the change in CO2 and uptake, we can calculate evasion.

# Distances are 17,30,80

#filt$evasion_14 <- ((filt$V1-(filt$V4+filt$Uptake))/ 100) *.0018*(1000000/44.01)*filt$stn3_Q/1000 #with 'uptake'

filt$evasion_14 <- ((filt$V1-filt$V4)/ 100) *.0018*(1000000/44.01)*filt$stn3_Q/1000 #without

Plot data

Histogram of Evasion Stations 1 -> 4

ggplot(filt, warning = FALSE, message=FALSE)+
  geom_histogram(aes(x = evasion_14),fill = "#767676",color="#4B9CD3",binwidth = 1)+
  geom_vline(xintercept = 0, linetype="dotted", color = "#4B9CD3", size=1.5)+
  xlim(-10,40)+
  ylim(0,2600)+
  labs(x = "Evasion", y="Count of values",title = "Station 1 -> 4")+
  theme_bw()
## Warning: Removed 2 rows containing missing values (geom_bar).

Evasion and time

ppt and Evasion

pptBar <- ggplot()+
  geom_bar(data = pptDaily,stat = "identity", aes(x = DateTime, y = ppt_mm))+
  geom_point(data = filt,aes(x = DateTime, y = evasion_14))


pptBar

##Evasion

dayPoly <- data.frame(x = c(as.POSIXct("06:15", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("18:21", format = "%H:%M"),as.POSIXct("06:15", format = "%H:%M")),y=c(-35,-35,40,40))

plot14 <- ggplot(filt)+
  geom_polygon(data = dayPoly,aes(x=x,y=y),fill="#f2e198")+
  geom_vline(aes(xintercept = as.POSIXct("06:15", format = "%H:%M")))+
  geom_vline(aes(xintercept = as.POSIXct("18:21", format = "%H:%M")))+
  geom_line(aes(x = time, y = evasion_14,group = day, col=stn1_Q),lwd=2)+
  scale_color_gradient(name = "Discharge (L/s)", trans = "log",
                        breaks = my_breaks, labels = my_breaks)+
  labs(x = "Time",y="Evasion", title = "Evasion station 1 -> 4")+
  theme_bw()

plot14

Evasion vs. Evasion

ggplot(filt)+
  geom_point(shape=21,aes(x = evasion_14, y=Flux_1, fill = "eosFD 1"), fill="#4B9CD3", alpha =.6,size=2)+
  geom_point(shape=21,aes(x = evasion_14, y=Flux_2, fill = "eosFD 2"), fill="#C7DA7D", alpha = .6,size=2)+
  geom_smooth(aes(x=evasion_14,y=Flux_1),method='lm',formula=y~x, color = "#4B9CD3")+
  geom_smooth(aes(x=evasion_14,y=Flux_2),method='lm',formula=y~x,color ="#C7DA7D")+
  xlim(0,36)+
  ylim(0,6)
## Warning: Removed 3394 rows containing non-finite values (stat_smooth).
## Warning: Removed 3379 rows containing non-finite values (stat_smooth).
## Warning: Removed 3394 rows containing missing values (geom_point).
## Warning: Removed 3379 rows containing missing values (geom_point).

Evasion based on Daily stream metabolism

Using the USGS streammetabolizer package, daily estimates of GPP and ER can be calculated at each of the DO sensor locations. This is done in a seperate script located in: ‘Ecuador/Analysis/Stream_Metabolism/metabolism.R’. The results are used to estimate daily evasion for the reach. The values expressed in the results of the streammetabolizer tool are grams of O2 / m^2 /d^-1. We need to convert this into uMols and then multiply by -1 to determine uMols CO2.

pred1 <- read.csv(here("Analysis/Stream_Metabolism/ModelOutputs/DO_1_Predictions.csv"))%>%
  select(date,GPP,ER)%>%
  mutate(Station = "1")
pred2 <- read.csv(here("Analysis/Stream_Metabolism/ModelOutputs/DO_2_Predictions.csv"))%>%
  select(date,GPP,ER)%>%
  mutate(Station = "2")
pred4 <- read.csv(here("Analysis/Stream_Metabolism/ModelOutputs/DO_4_Predictions.csv"))%>%
  select(date,GPP,ER)%>%
  mutate(Station = "4")

metab <- rbind(pred1,pred2,pred4)%>%
  mutate(DateTime = as.POSIXct(date))

metab%>%
  plot_ly(x = ~DateTime, y = ~GPP,hoverinfo = "text",
        text = ~paste("Station: ", Station,
        "<br> Date: ", DateTime,
        "<br> ER: ", round(ER,2)))%>%
  add_markers(color = ~Station)%>%
  layout(
    title = "Gross Primary Productivity (DO)",
    xaxis = list(title = "Date"),
    yaxis = list(title = "GPP (g/m^2/d^-1)")
  )
## Warning: Ignoring 42 observations

Estimate daily ER and GPP in uMols of CO2

1 mole of CO2 = 31.9988 grams. We average the GPP and ER values from the three stations to estimate GPP and ER for the entire reach.

metabCO2 <- metab%>%
  mutate(GPP_CO2uMol = GPP/31.9988*-1000000,
         ER_CO2uMol = ER/31.9988*-1000000)%>%
  group_by(DateTime)%>%
  mutate(AvgGPP_CO2 = mean(GPP_CO2uMol),
         AvgER_CO2 = mean(ER_CO2uMol))

metabCO2%>%
  plot_ly(x = ~DateTime, y = ~GPP_CO2uMol,hoverinfo = "text",
        text = ~paste("Station: ", Station,
        "<br> Date: ", DateTime,
        "<br> ER: ", round(ER_CO2uMol,2)))%>%
  add_markers(color = ~Station)%>%
  layout(
    title = "Gross Primary Productivity (uMols CO2)",
    xaxis = list(title = "Date"),
    yaxis = list(title = "GPP (uMols/m^2/d^-1)")
  )
## Warning: Ignoring 42 observations
metabCO2%>%
  plot_ly(x = ~DateTime, y = ~ER_CO2uMol,hoverinfo = "text",
        text = ~paste("Station: ", Station,
        "<br> Date: ", DateTime,
        "<br> ER: ", round(GPP_CO2uMol,2)))%>%
  add_markers(color = ~Station)%>%
  layout(
    title = "Ecosystem Respiration (uMols CO2)",
    xaxis = list(title = "Date"),
    yaxis = list(title = "ER (uMols/m^2/d^-1)")
  )
## Warning: Ignoring 42 observations

Take our evasion estimates and sum them by day for comparison

evasionDaily <- filt%>%
  group_by(day)%>%
  mutate(DailyFlux = mean(evasion_14*86400))%>%
  select(day,DailyFlux)%>%
  distinct()%>%
  rename("date"="day")

Merge ‘Daily Evasion’ with stream metabolism to get metabolism adjusted evasion

metabAvg <- metabCO2%>%
  select(date,AvgGPP_CO2,AvgER_CO2)%>%
  distinct()
## Adding missing grouping variables: `DateTime`
adjustEvasion <- metabAvg%>%
  left_join(evasionDaily,by="date")%>%
  mutate(AdjEvasion = DailyFlux+AvgGPP_CO2+AvgER_CO2)%>%
  select(-date)
## Warning: Column `date` joining factors with different levels, coercing to
## character vector
adjustEvasion%>%
  plot_ly(x = ~DateTime, y = ~AdjEvasion,hoverinfo = "text",
        text = ~paste("Evasion: ",round(AdjEvasion,1),"uMols CO<sub>2</sub> m<sup>2</sup>d<sup>-1</sup>",
                      "<br> GPP: ", round(AvgGPP_CO2,1),"uMols CO<sub>2</sub> m<sup>2</sup>d<sup>-1</sup>",
        "<br> ER: ", round(AvgER_CO2,1),"uMols CO<sub>2</sub> m<sup>2</sup>d<sup>-1</sup>"))%>%
  add_markers(color = "Evasion", size=50)%>%
  add_markers(y = ~AvgER_CO2, color = "Respiration", size = 50)%>%
  add_markers(y = ~abs(AvgGPP_CO2), color = "Photosynthesis", size=50)%>%
  layout(
    title = "Daily Evasion (Entire Reach)",
    xaxis = list(title = "Date"),
    yaxis = list(title = "uMols m<sup>2</sup>d<sup>-1</sup>", type = "log")
  )
## Warning: Ignoring 15 observations
## Warning: Ignoring 14 observations

## Warning: Ignoring 14 observations